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We present a Monte Carlo approach to incorporating the effect of thermal fluctuations 
in field theories of polymeric fluids. This method is applied to a field-theoretic model of 
a ternary blend of AB diblock copolymers with A and B homopolymers. We find a shift 
in the line of order-disorder transitions from their mean-field values, as well as strong 
signatures of the existence of a bicontinuous microemulsion phase in the vicinity of the 
mean- field Lifshitz critical point. This is in qualitative agreement with a recent series of 
experiments conducted with various three-dimensional realizations of this model system. 
Further, we also compare our results and the performance of the presently proposed 
simulation method to that of an alternative method involving the integration of complex 
Langevin dynamical equations. 



I. INTRODUCTION 



A straightforward approach to improving the 
properties of polymeric materials, such as their 
stiffness, conductivity, toughness, is to blend two 
homopolymer species (A and B) into a single 
melt 1,3 ' 5,4 . However, in most cases, weak segre- 
gation tendencies between A and B monomers 6 
usually cause such a direct blend to macroscop- 
ically phase separate below a critical tempera- 
ture, or, equivalently, above a critical strength of 
interaction 2 . Such phase separation often leads to 
poor quality materials with non-reproducible mor- 
phologies and weak interfaces. To circumvent this 
problem, one typically introduces compatibilizers 
into the system, most commonly multicomponent 
polymers like block, graft, or random copolymers 7 . 
Such compatibilizers preferentially segregate to- 
wards the interfaces between the bulk A and B do- 
mains, thereby lowering interfacial tensions, stabi- 
lizing the formation of complex morphologies 8-11 , 
and strengthening interfaces 12-15 . The compati- 
bilizers so added are usually the most expensive 
component of a blend, and consequently, theoret- 
ical studies of the phase behavior of such blends 
possess important ramifications in rendering such 
industrial applications more economically feasible. 

Here we examine the phase behavior of a ternary 
blend of symmetric A and B homopolymers with 
an added AB block copolymer as a compatiblizer. 
Specifically, we consider symmetric diblock copoly- 
mers wherein the volume fractions of the two dif- 
ferent blocks in the copolymer are identical. More- 
over, the ratio of the homopolymer and copolymer 
lengths is denoted as a and is chosen as a = 0.2. 



This is largely motivated by a recent series of 
experiments by Bates et al. 16-19 on various re- 
alizations of the ternary AB+A+B system with 
a rs 0.2. A combination of neutron scattering, 
dynamical mechanical spectroscopy, and transmis- 
sion electron microscopy (TEM) was used to ex- 
amine the phase behavior. An example of an ex- 
perimental phase diagram is shown in Fig. 1. It 
will be discussed in detail further below. 

Theoretically, one of the most established ap- 
proaches to studying such polymer blends is the 
self-consistent field theory (SCFT), first intro- 
duced by Edwards 20 and Helfand 21 and later 
used, among others, by Scheutjens and Fleer 22 ' 23 , 
Noolandi et al. 24 , and Matsen et al. 10 . In imple- 
menting SCFT 25 , one typically chooses (a) a model 
for the polymers and (b) a model for the interac- 
tions. Most researches adopt the Gaussian chain 
model, wherein the polymers are represented as 
continuous paths in space 26 . It models polymers 
as being perfectly flexible, i.e., there exists no free 
energy penalty for bending. Instead, it features a 
penalty for stretching. (In contexts where the lo- 
cal rigidity or orientation effects are important, the 
wormlike chain model 27 serves as a popular alter- 
native.) Such a molecular model for the polymers 
is typically supplemented by a model for the in- 
teractions between unlike polymers, which is usu- 
ally chosen as an incompressibility constraint plus 
a generalization of the Flory-Huggins local contact 
interaction term, whose strength is governed by the 
product of the segregation parameter, %j and the 
polymerization index, N . 
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FIG. 1. Experimental phase diagram 19 of a 
PDMS-PEE/PDMS/PEE blend with a « 0.2. Repro- 
duced by permission of The Royal Society of Chem- 
istry. 



We have implemented SCFT for the above 
model of polymers, and investigated the phase be- 
havior of the ternary A+B+AB blend. Fig. 2 
displays the resulting mean-field phase diagram 
in the coordinates of homopolymer volume frac- 
tion 4>h and segregation strength %A '. Note that 
the phase diagram shown is a slice along the "iso- 
plcth" where the two homopolymers are blended in 
equal proportions. It can be observed that on the 
copolymer-rich side of the diagram shown in Fig. 
2, a line of order-disorder transitions (ODT) sep- 
arates the disordered region (which occurs at low 
interaction strengths, \N) from a periodically or- 
dered (i.e., lamellar) phase at higher %iV. For pure 
copolymers (<pu = 0), this transition is predicted 
within mean-field theory to occur at \N = 10.495, 
and the line of ODTs, often referred to as Leibler 
line, represents second-order transitions. On the 
homopolymcr-rich side of the diagram, we find a 
line of second-order transitions from a disordered 
phase at low xN to a region of two coexisting ho- 
mogeneous liquid phases at high \N ■ For the pure 
homopolymer system (<j>H = 1), the upper criti- 
cal consolute point is found by mean-field Flory- 
Huggins theory to be at a\N = 2. This line of 
continuous transitions is referred to as Scott line 28 . 
The point where the Leibler and Scott lines meet 
is found, again within mean-field theory, to be an 
isotropic Lifshitz critical point (LP) 29 , which in the 
case of a = 1 becomes a Lifshitz tricritical point 
(LTP) according to Broseta and Fredrickson 30 , and 
is predicted to occur at a total homopolymer frac- 
tion of <f>L = 1/(1 + 2cv 2 ) and an incompatibil- 
ity of {xN) L = 2(1 + 2a 2 ) /a. The separation 
into microphases (i.e., lamellae) observed along the 
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FIG. 2. Mean-field phase diagram for a ternary 
AB+A+B blend with a = 0.2, as obtained from a 
Fourier space implementation of SCFT. 2<j) denotes a 
region of two-phase coexistence between an A-rich and 
a B-rich phase, 30 one of three-phase coexistence be- 
tween an A-rich, a B-rich, and a lamellar phase. See 
text for more explanation. 

lamellar periodicity until it finally diverges at the 
LP, giving rise to macrophase separation along the 
Scott line. Below the LP in the ordered regime, 
the system undergoes a first-order transition from 
the lamellar (L) to the two-phase (A+B) region 
along the axis of homopolymer concentration, <pu- 
Contrary to a previous publication 18 in which this 
transition was displayed as continuous, a care- 
ful reexamination of the mean-field theory reveals 
that it is indeed first-order with a three-phase 
(L+A+B) coexistence region reaching all the way 
up to the LP for any a. This is further corrobo- 
rated by recent results of Naughton and Matsen 31 . 

The mean-field phase diagram displayed in Fig. 
2 differs substantially from the experimental phase 
diagram, shown in Fig. 1. In the experiments, the 
LP was destroyed by fluctuations, giving rise to 
a bicontinuous microemulsion phase in its vicin- 
ity. At lower temperature T, the L and A+B 
phases were conjectured to be separated by a nar- 
row channel of this microemulsion phase stretching 
down from the Lifshitz region, while the Leibler 
and Scott lines were shifted from their mean-field 
locations to lower temperatures (i.e., higher %./V, 
which is proportional to 1/T). Experimental re- 
sults are consistent with the picture that as the 
lamellar periodicity increases along the 0# axis, 
the persistence length of composition fluctuations 
along the microdomain boundaries decreases un- 
til at some point the two become comparable in 
size and the lamellae begin to rupture forming a 
bicontinuous structure. 

The experimental results discussed above, while 
in contradiction to the mean-field SCFT calcu- 



lations, are not entirely surprising when viewed 
in the context of theoretical studies pertaining to 
the effect of fluctuations on the mean-field phase 
diagrams. Indeed, even in the context of the 
pure diblock copolymer (i.e. <pn = 0), Leibler 32 
speculated that the order-disorder transition (pre- 
dicted to be second order within the mean-field 
theory) would become a weakly first-order transi- 
tion when the effects of fluctuations are accounted 
for. Such a reasoning was confirmed quantitatively 
by Fredrickson and Hclfand 33 , who extended an 
earlier analysis by Brazovskii 34 to show that fluc- 
tuations change the order of the transition and also 
shift the transition to higher incompatibilities, X-/V. 
Further, they quantified such shifts in terms of a 
parameter that depends on the polymer density 
(or, cquivalcntly in the context of a polymer melt, 
their molecular mass N). This parameter, denoted 
generically as C in this text, acts as a Ginzburg 
parameter 35 , and in the limit C — > oo, the mean- 
field solution is recovered. 

The present study is motivated by the hypoth- 
esis that the fluctuation effects not accounted for 
within SCFT are responsible for the discrepancies 
noted between the experimental results and the 
theoretical predictions for the ternary-blend phase 
behavior. As a partial confirmation of our hypoth- 
esis, the experimental phase diagrams (Fig. 1) of 
Bates and co-workers deviate most clearly from 
the mean-field diagram for the PEE/PDMS/PEE- 
PDMS system, in situations where the molecular 
weights were much lower than in the other systems 
(i.e., where the C parameter is small). This arti- 
cle presents a theoretical examination of the effect 
of thermal fluctuations upon the phase behavior 
of ternary blends. Our primary goal is to under- 
stand and quantify the shift in both the Leibler and 
Scott lines and to examine the formation of the mi- 
croemulsion phases in such systems. Bicontinuous 
structures such as the ones observed in the above 
experiments are particularly interesting from an 
application point of view in that they are able to 
impart many useful properties on polymer alloys. 
For example, if one blend component is conducting 
or gas-permeable, these properties will be passed 
on to the entire alloy. Also, mechanical proper- 
ties such as the strain at break and the toughness 
index have been observed to exhibit maxima with 
co-continuous structures 36 . 

Previous simulation studies of fluctuation ef- 
fects in ternary diblock copolymer / homopoly- 
mer blends 37 ' 38 ' 41 ' 39 ' 40 , have mostly focused on the 
a = 1 case. Within a lattice polymer model, 
Muller and Schick 37 have investigated the phase di- 
agram of a corresponding symmetrical A+B+AB 
blend. According to their results, the disordered 
phase extends to higher %iV than predicted by the 



mean field theory, and the LTP is replaced by a 
regular tricritical point. Due to finite size effects, 
the transition to the lamellar phase could not be 
localized quantitatively. 

Particle-based models have the disadvantage 
that one is restricted to studying relatively short 
polymers. This motivates the use of field theo- 
ries in simulations. Kiclhorn and Muthukumar 41 
have studied the phase separation dynamics in the 
vicinity of the LTP within a Ginzburg-Landau free 
energy functional, which they derived from the Ed- 
wards model 20 within the random phase approxi- 
mation (RPA) , in the spirit of the theories for pure 
diblock copolymer systems mentioned earlier 32,33 . 
Naively, one should expect that such a description 
is adequate in the weak segregation limit, where 
the correlation length for composition fluctuations 
is large. However, Kudlay and Stepanow 42 have re- 
cently questioned this assumption. They pointed 
out that the fluctuations in RPA-derived Hamilto- 
nians have significant short wavelength contribu- 
tions which cannot be eliminated by rcnormaliza- 
tion and may lead to unphysical predictions. 

In this work, we propose an alternative field- 
theoretic approach which rests directly on the Ed- 
wards model, without resorting to the additional 
RPA approximation. The partition function is 
evaluated with a Monte Carlo algorithm. A re- 
lated algorithm, the complex Langcvin method, 
has been presented recently by two of us 43-45 . 
Here we study two-dimensional ternary blends at 
a = 0.2 with both methods and compare the re- 
sults. 

The remainder of this article is organized as fol- 
lows. In section II, we outline the field-theoretic 
formulation of our model. The simulation method 
is presented in section III. In section IV, we 
present and discuss our results before we conclude 
in section V with a summary and an outlook on 
future work. 



II. THE FIELD-THEORETIC MODEL 

In this section, we present a field-theoretic model 
for the system under consideration. To this effect, 
we consider a mixture of ua homopolymers of type 
A, riB homopolymers of type B, and uab symmet- 
ric block copolymers in a volume V. The polymer- 
ization index of the copolymer is denoted N, and 
the corresponding quantities for the homopolymers 
are denoted Na = Nb = ceN. All chains are as- 
sumed to be monodisperse. We consider the case 
of symmetric copolymers, where the fraction of A 
monomers in the copolymer is / = 1/2. We re- 



strict our attention to the concentration isoplcth, 
where the homopolymers have equal volume frac- 
tions <j>HA — <I>hb — 4>h/2. The monomeric vol- 
umes of both A and B segments are assumed to 
be identically equal to 1/po- In this notation, the 
incompressibility constraint requires 

4>hPqV paV 
n A = n B = - Ar ; n AB = (1 - 4>h)— rp, (1) 



2aN 



N 



where V denotes the volume of the system. In the 
Gaussian chain model, the segments are assumed 
to be perfectly flexible and are assigned a common 
statistical segment length b. The polymers are rep- 
resented as continuous space curves, R}(s), where 
i denotes the polymer species (A, B, or AB), j 
the different polymers of a component, and s rep- 
resents the chain length variable measured along 
the chain backbone such that < s < uj with 
vab = 1 and v A = v B = ot. Conformations of non- 
interacting polymers are given a Gaussian statisti- 
cal weight, cxp(— Hq), with a harmonic stretching 
energy given by (units of k B T), 



Furthermore, the incompressibility constraint is 
implemented by the insertion of a Dirac delta func- 
tional 46 , and we obtain the canonical partition 
function of the blend, Zc- 

Z c oc f nil^ R }( s ) cxp(-H - Hj), 

J » 3 

(7) 

where the J DR denote functional integrations 
over chain conformations. Representing the delta 
functional in exponential form and employing a 
Hubbard-Stratonovich transformation on the rh 2 
term in (4), we obtain 



Z c cx f DW- f DW+cxp[-H c (W-,W + ) 

J oo J zoo 



with 

H C (W-,W+) = C 
-V{l-<f> H )\nQ AB 



(8) 



1 



drWl 



drW+ 



V4>h 
2a 



(\nQ A + \nQ B ) (9) 



dR)(s) 



ds 



(2) 



where b A = b B = b AB = b is the statistical seg- 
ment length. Interactions between A and B seg- 
ments are modeled as local contact interactions by 
a pseudopotential with a Flory-Huggins quadratic 
form in the microscopic volume fractions, 



H, 



[R] = xpo J 



dr cj> A (r)(f) B (r) 



(3) 



which with the substitutions rh — <j) B — 4> A , 
<Pa+ <I>b transforms into 



H,\R 



dr ((j) 2 -rh"). (4) 



X denotes the Flory-Huggins parameter, and the 
microscopic density operators <f> A and <f) B are de- 
fined by 

A (r) = -£ / ds5(r--Rf B (s)) 
Po £-f Jo 



m nA r a 
-J2 / dsS(r-Rf(s)) (5) 
Po ^ Jo 



JV nAB f 1 
Po j=1 Jf 



N 



»B pa 

]T / dsS(r-Rf(s)). (6) 



P0 J 



r - Po R d 



(10) 



Here and throughout this article, all lengths are ex- 
pressed in units of the unperturbed radius of gyra- 
tion, R go = b(N/(2dj) 1 / 2 , where d is the space di- 
mension. The parameter C in the above equations, 
which occurs as a global prefactor to He, acts as a 
Ginzburg parameter such that in the limit C — > oo 
the partition function (8) is dominated by its sad- 
dle point and the mean-field solution becomes ex- 
act. In (9), i = V^T, and Q A ,Q B , and Q AB de- 
note respectively the single-chain partition func- 
tions for the A, B, and AB chains, in the potential 
fields VF_(r) and W + (r). Note that W- is conju- 
gate to the difference in A and B densities, rh, and 
W + to the total density, <j>. Moreover, W_ is real, 
whereas W+ is imaginary, thereby rendering He 
complex. 

The single-chain partition functions can be ex- 
pressed in terms of the Feynman-Kac formulae 21 
as: 



Qi = / dr qi{r,Ui), 



(11) 



where the propagators qi satisfy the diffusion equa- 
tions 

d 

—q l (r,s)=Aq l (r,s)-q i U l ; ft(r,0) = l (12) 



U A = W+-W-, 0<s<a 
U B = W+ + W-, 0<s<a 

'U A , 0<s<f 
U B , f<s<l 



(13) 



Vab = 



An analogous diffusion equation applies to the con- 
jugate propagators, q\, which propagate from the 
opposite end of a polymer. Due to their symme- 
try, the propagators of the homopolymers are their 
own conjugates. Hence, we can calculate density 
operators, (f> A and (f> B , from 

<M r ) = 77— / ds lAB{r, s) q AB (r, 1-s) 



!AB Jo 



V6 H f a 

+ j\ \ ds q A (r,s) q A (r,a- s), (14) 

V (1 (j) ) f ^ 

fe(r) = ^ — / ds q AB (r, s) q AB {r, 1-s) 

WAB J f 



2aQ B 



/ ds q B (r, s) q B (r, a - s). (15) 
Jo 



These d shown on these snapshots are not to be 
confused with real densityensity operators depend 
on W- and W+ and arc in general complex. How- 
ever, their ensemble averages yield real densities 
which correspond to the experimentally measur- 
able quantities: <j> A ,B = (4> AtB }(= {4>a,b})- 

The above derivation assumed a canonical en- 
semble for the blend, wherein the average compo- 
sitions of the different compositions are fixed. For 
future reference, we also list a similar set of formu- 
lae obtained in the grand-canonical ensemble: 

H GC (W-,W+) = C J dr W 2 _ - j dr W+ 



-Qab — zQa — zQ B 



(16) 



where z = exp((A/j,/k B T),a,nd A/i = [i A — jj, AB = 
f-B — Hab is the difference in the chemical poten- 
tials of homopolymers and copolymers. The den- 
sity operators are given as: 

4> A (r) := j dsq AB (r, s)q AB (r, 1-s) 
Jo 

+z ds q A (r,s)q A (r,a - s), (17) 
Jo 

4> B {r) := ds q AB (r, s)q AB (r, 1 - s) 

+z ds q B (r, s)q B (r, a - s). (18) 
Jo 

III. FIELD-THEORETIC SIMULATIONS 
A. General considerations 



Field-theoretic models like the one introduced 
above have been used in many earlier researches, 



and the resulting functional integrals of the parti- 
tion function have been analyzed using various ap- 
proximate methods. In self-consistent mean-field 
theory (SCFT) 10 ' 25 ' 47 , the integrand of (8) is ap- 
proximated by saddle points (Wf , W*) such that 



5H 



5W+ 



0, 



5H 



0. (19) 



A homogeneous saddle point then corresponds 
to the disordered phase, whereas every ordered 
phase has a unique inhomogeneous saddle point. 
By examining (8), it is further apparent that C 
plays the role of a Ginzburg parameter in the 
sense that the mean-field approximation, i.e. Z ~ 
exp(— H[W+, W*]), becomes exact in the formal 
limit of C — > 00, corresponding to N — > 00 for 
d > 2. Implementations of SCFT 10 can thus 
be viewed as numerical strategies for computing 
the various saddle points. However, thermal fluc- 
tuation effects in the vicinity of phase bound- 
aries invalidate the mean-field results in these re- 
gions, and strategies for evaluating (8) must be 
found. Analytical methods are quite limited, es- 
pecially for systems dominated by inhomogeneous 
saddle points, and often involve weak segrega- 
tion expansions . In contrast, earlier publica- 
tions of two of us 43 ~ 45 have detailed a new ap- 
proach termed field-theoretic simulations, which is 
a methodology for numerically sampling functional 
integrals in field-theoretic models of polymer solu- 
tions and melts. The underlying idea is to nu- 
merically sample Z (independently of the type of 
ensemble employed) by generating a Markov chain 
of W- , W + states with a stationary distribution 
proportional to exp(— if [VIA, W + \). Averages of 
physical quantities f[W- , W+] can then be approx- 
imated by "time" averages using the states of the 
Markov chain, 

M 

(/[HA , W+]) « 1/M £ f[W-j , W+j], (20) 
j=i 

where j labels the states of the Markov chain. 
True equality of ensemble and time averages is 
established as usual for ergodic systems in the 
limit M — » 00. Such a numerical method is non- 
perturbative in nature, allowing for a more com- 
plete account of fluctuation effects. 

In this article, we undertake field-theoretic sim- 
ulations to address the effect of fluctuations upon 
the mean-field phase diagram. Note that to nu- 
merically generate the Markov chain of HA , W + 
states, one needs to account for the fact that 
if[W_,W+] ((9) and (16)) in general possesses 
both real and imaginary parts, implying a non- 
positive-definite statistical weight exp(-H). In 



our previous researches we implemented a gen- 
eral "complex Langevin dynamics" technique that 
had been proposed earlier to deal with such non- 
positive-definite statistical weights 48 ' 49 . In view of 
the fact that this method has been discussed in 
detail in our previous publications, we restrict our 
discussion to a brief outline of the technique which 
is presented in the Appendix. It is to be noted 
that outstanding issues still remain as regards the 
theoretical foundations of this technique 48 . While 
in our earlier studies we encountered no problems 
with the convergence or the uniqueness of the so- 
lutions, in the present study this method did en- 
counter numerical difficulties, especially in the re- 
gions close to the Lifshitz point. These difficulties 
are associated with strong fluctuations and com- 
plex phase oscillations and will be described more 
fully in a subsequent section. 

In an effort to probe and overcome the difficul- 
ties encountered in preliminary complex Langevin 
simulations of the present ternary blend and quan- 
tify the effect of fluctuations on the phase diagram 
of this system, we present an alternative, albeit 
equally rigorous technique for implementing field- 
theoretic simulations, which we refer to as the field- 
theoretic Monte Carlo method. In the discussion of 
our results, we compare explicitly the results ob- 
tained by the latter method with that obtained 
through the complex Langevin method and find 
good agreement in the range of their validities. 
It is pertinent to also note that there also exist 
real Langevin methods in the literature that focus 
specifically on the W- fluctuations 50 . 



B. The field-theoretic Monte Carlo method 

Our approach to sampling (8) is based on the 
Monte Carlo method. We recall that W_ is real 
and W + imaginary, and that the fluctuations in 
W + are conjugate to the total density. Although 
the incompressibility constraint implemented in 
our model effectively only constrains <j> = cpA + 
4>b = (4>a + 4>b), one might guess that the effect of 
the W+ fluctuations that impose this constraint is 
rather small compared to the composition fluctua- 
tions governed by W-. Therefore, we shall tackle 
W- and W + in two distinct steps, after which we 
will be able to argue that the effect of the W- 
fluctuations is indeed the dominant one. 

To simulate the fluctuations in W-, we pick a 
starting density configuration and use the SCFT 
mean-field equations 25 to calculate a seed W- and 
W+. Next, we proceed according to the following 
scheme: 



(a) Make a global move in W-, i.e., tentatively 
add a random number G [— 1,+1] times a step 
width parameter to every W-(r), and denote by 
W° the old 

(b) Calculate a self-consistent W + that is a par- 
tial saddle point with the tentative new W- fixed; 
to do this, we need to iteratively solve the equation 

Mr; [W+, W-]) + Mr; [W + , W-]) = 1 (21) 

for W+(r). The initial guess for W+ in this iter- 
ation, W®, is taken to be the old W+ from the 
last Monte Carlo cycle. The iteration, which em- 
ploys a two-step Anderson mixing scheme 51 , is ter- 
minated once AH = H[W-, W+ } - H[W°,Wl\ is 
determined to within 0.001, or a few parts in 10 4 . 
Every step in this iteration requires the solution 
of the diffusion equation (12) on the entire simu- 
lation lattice, and this is where practically all the 
computing time is spent. The number of iterations 
to find W+ in most cases is roughly around 10. 

(c) Accept or reject the move according to a 
standard Metropolis criterion. In other words, we 
always accept the move if the resulting difference 
in H, i.e., AH, is negative, and otherwise, we ac- 
cept with a probability exp(— AH). 

(d) Go back to (a) unless the maximum number 
of Monte Carlo steps has been reached. 

Next, we turn to discuss the sampling of W + 
fluctuations. As a result from SCFT, we know that 
the saddle point of W+ as well as the partial 
saddle point W+[VF_] are purely real, despite the 
fact that the integration path of W+ in (8) is along 
the imaginary axis. Neglecting surface terms and 
because we do not cross any singularities, we can 
therefore deform the integration path and repre- 
sent W + as 

W+(r) = W+[W-]{v) +iCu+(r), (22) 

where cD+(r) is real. Furthermore, we split the ar- 
gument of the integral in the partition function 
into a real part and a complex reweighting factor, 
leading to the following expression for the canoni- 
cal ensemble: 

Z c oc J VW. J VCo+ cxp J c , (23) 



drWl 



drW+ 



(24) 



Ic 



= exp (c \i J 



AvCj i 



+ t^^Im(lnQ i [Wl,a> + ])l). (25) 



TjR 



where V% is the volume occupied by species i. In 
the grand canonical ensemble the corresponding 
expression is: 

Z GC cx J VW_J V£o + exp (-H§ c ) ■ I GC , (26) 

_ c[±JdrWl-JdrW + 

-^z i Re(Q i [W-,£j+]))], (27) 

i 

Igc = cxp (c \i J drw+ 

+ i S ^z l Im(Qi[W_,w + ])]), 



where zab — 1 and za = zb = z. 



We can thus in principle simulate both W- and 
W + fluctuations, e.g., by alternating a move in 
W- with a number of moves in uj + and replacing 
AHc.gc with AHq GC in the Metropolis criterion. 
Each configuration then needs to be weighted with 
the factor Ic,GC for the purpose of computing av- 
erages. We note that the propagators and the field 
term in the diffusion equation (12) are now com- 
plex, but the averages of all physical quantities will 
approach real values in equilibrium. 

At this point, a remark concerning our numerical 
methods is due: the single most time-consuming 
part in doing our simulations is the solution of 
the diffusion equation (12). Initially, we used ei- 
ther the Crank-Nicholson (CN) or DuFort-Frankel 
(DF) finite differencing schemes 52,53 to effectuate 
this task. As it turned out, however, a pseudo- 
spectral method that had been used successfully 
in the past to solve the mathematically very sim- 
ilar Schrodingcr equation and that had recently 
been applied to SCFT 54 ' 55 was far better suited 
for our needs. Typically, it allowed for an order of 
magnitude higher values of contour step ds than 
cither CN or DF to achieve the same accuracy. 
Implementing the scheme with highly optimized 
Fast Fourier Transform routines like the publicly 
available FFTW library 56 , the computing time per 
contour step is only roughly double that needed for 
DF and half that needed for CN on a single pro- 
cessor. We therefore regard it as the method of 
choice in the context of polymer field theory. For 
the larger system sizes to be discussed in the next 
section, a parallel version of the code was used. 



C. Analysis of W+ fluctuations 



the complex reweighting factor Ic or Igc an d (b) 
through the difference of AH R from AH. 

We consider first Igc and to this end split the 
complex diffusion equation into two real equations 
for the real and imaginary parts of the propagator 
for polymers of type i, denoted gf and qf , respec- 
tively. Note the symmetry properties of this cou- 
pled set of equations with respect to a sign change 
in u> + : for W + = W + ± iui + we obtain 



d_ R _ 

9 s % 



Agf - U iq f- ± u+ql 
± M = ± Aq{ T U t ql - u, + q* 



(29) 



(28) We see that qf- remains unchanged whereas q\ goes 
to —q\. Thus, for the corresponding single-chain 
partition function Qi = Qf + iQl , Qf is symmet- 
ric in Q + , and Q\ antisymmetric. Furthermore, 
the exponent of the reweighting factor Igc is anti- 
symmetric in (2)+ (cf. (28)). Since in the partition 
function (26) we integrate over all possible config- 
urations of iZ>_|_, the imaginary parts of Igc cancel 
out for each antisymmetric pair, and consequently, 
we effectively only need to apply the real part of 
Igc'- 



GC 



cos C 



J drw+ + '%2z i Q![W-,u> + ] 



(30) 

These symmetries reveal that in a functional 
expansion of Qi in uj + around the partial saddle 
point, the even-order terms coincide with 

the expansion of Qf , and the odd-order terms with 
that of Ql . Moreover, 



SQj = 6Q! 
SW+ 5ib + 



and 



SQj 
5W+ 



-05) 



(31) 



(32) 



by virtue of the incompressibility constraint. 
Hence, 



+ 0(u>l), (33) 



and finally 

Igc = cos (0(u% )) = 1 - 0{Co\), (34) 
Qf = Q t [WA+0(u 2 + ). (35) 



Evidently, the W + fluctuations exert their influ- 
ence on the simulation in two ways: (a) through 



Based on these analytical results, we should ex- 
pect Igc to deviate only very slightly from unity 
so that it can be neglected in most practical cases. 



Hence, the impact of the W + fluctuations would 
materialize solely as a (real) 0(Q+) contribution to 
Hqc- The above derivation of (34) and (35) was 
for the grand canonical ensemble but holds true if 
one is to derive analogous relations in the canon- 
ical ensemble, as well, because the leading term 
of hi Qi is linear in lu + , as in the grand canon- 
ical ensemble. It should however be noted that 
in the canonical ensemble, the gauge invariance of 
the theory to a constant shift of the <D+ field al- 
lows the introduction of a spurious phase factor 
e 1 ^ that preserves He (and hence \Qi\ and H R ). 
This if) field will freely move in canonical simula- 
tions unless constrained, e.g., by demanding that 
J dnl> + (r) = 0. In contrast, in the grand canonical 
ensemble, Q R and hence H R are to lowest approx- 
imation quadratic functions around Co + = 0. As a 
result the shape of the H R (W+) energy landscape 
is much steeper in the grand canonical ensemble 
than in the canonical, and no constraint is needed. 

While the above considerations can only give us 
an indication of the insignificance of the W+ fluc- 
tuations, we have also done test simulations to get 
more conclusive estimates. We therefore choose 
the grand canonical ensemble for testing our simu- 
lation technique of W+ fluctuations. Before start- 
ing a run at z = 2, %A = 12, a = 0.2, and C = 50, 
we fully equilibrated a simulation of W- fluctua- 
tions only with the same parameters, which started 
from a disordered configuration and ran a million 
Monte Carlo steps. This set of parameters is very 
close to the (fluctuation-corrected) order-disorder 
transition in the phase diagram, on the disordered 
side. The average homopolymer fraction, (j>jj, is 
0.65. On a 32x32 lattice with periodic boundary 
conditions, we used dx = 0.78125. We then took 
this equilibrated configuration and switched on the 
W + fluctuations, again making ten (global) moves 
in W + after each move in W- so as to roughly 
spend equal amounts of computing time for the 
two types of moves. Note in Fig. 3 how Hoc 
quickly increases and then plateaus off, indicating 
that the W + fluctuations have become saturated. 

For the same run, we recorded the quantity 
AH* = AH R - AH (cf. Eqs. (24) and (27)), 
which is dis played in t he inset of Fig. 3. One can 
interpret \/ ((AH*) 2 } as a measure for the devi- 
ation of the simulation from a reference case in 
which only W- fluctuations are simulated. This 
simulation was done with a chain-contour dis- 
cretization of ds — 0.01, which allows for the 
calculation of AH to within 0.01. Typically, in 
simulations involving the W- field alone, we used 
ds = 0.05 in order to save computing time. The 
discretization error of A H in this case is increased 
to < 0.08. Considering that typical values of AH 
are on the order of ±3, this appears like a very 
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FIG. 3. Ro(H) in a simulation of W- and W+ vs. 
number of Monte Carlo steps in W- for C = 50, 
xN = 12, z = 2 (grand canonical ensemble), started 
from a configuration equilibrated in W-. Each W- 
step alternates with ten steps in w+ . Inset shows AH* 
(difference of AH and a reference value at w+ fluctua- 
tions switched off) in the same simulation. 

moderate numerical error. At the same time, we 
realize that the effect of the fluctuations in W+ 
expressed through AH* is well below our stan- 
dard numerical accuracy: y/ ((AH*) 2 ) w 0.022. 
Furthermore, the average of AH* very closely ap- 
proximates zero, so there is no drift of A H as a 
result. We therefore conclude that contributions 
of the W+ fluctuations to H R cannot be distin- 
guished as such from numerical inaccuracies. 

Next, we examine the reweighting factor, Iqc-, 
in the particular simulation discussed above. Nu- 
merically, we must evaluate: 



I§ c = cos [C 



dru + + ZiQl ] = cos(A) 

(36) 

The two terms in the argument, A, of the co- 
sine in (36) have opposite signs and largely can- 
cel out each other, as was demonstrated theoret- 
ically above. However, in practice, the numerical 
application of Iqc to the simulation is impeded 
greatly by the fact that the two contributions to 
A are extensive, so the fluctuations in A are am- 
plified. Iqc is a non-positive-definite weighting 
factor, echoing the well-known "sign problem" en- 
countered in fermionic field theories in elementary 
particle physics 57 . Fig. 4 displays part of a time se- 
ries of A in our above test run. It performs fast os- 
cillations around a very small absolute mean value 
which within numerical accuracy closely approxi- 
mates zero when averaged over the entire run, al- 
though it is expected from theory to be finite. This 
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FIG. 4. The argument, A, of the cosine in the 
reweighting factor in the same simulation as in Fig. 3. 
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FIG. 5. Autocorrelation function of W- in the same 
simulation as in Fig. 3. 



is the essence of our sign problem: the statistical 
deviation from the mean completely dominates A. 
At this point, we seem to have arrived at a dead 
end, having to concede that we cannot numerically 
implement the reweighting factor. 

However, we still may get a sense of the impact 
of the uj + fluctuations from examining pairwise 
correlation functions of u> + , W-, and Iqc- The 
correlation function is defined as 



C[A,B;t] 



(A(t')B(t' + t))-(A)(B) 



2 Jl 



(37) 



A U B 



where a denotes standard deviation, and t and t' 
are numbers of moves in W- . If A or B is a field, 
we average over all lattice points r. Fig. 5 shows 
C[W-,W-}, and in Fig. 6, C[I,I], C[Q+,u)+], 
C[I,W-}, and C[u)+, W-] are displayed. On com- 
paring the various correlation times, an important 
qualitative difference is observed: while the au- 
tocorrelation time of VF_ is on the order of 10 5 
Monte Carlo steps in W- (subsequently denoted 
MCS), those of cD+ and / are only on the order of a 
few hundred MCS, i.e., several orders of magnitude 
smaller 58 . Furthermore, the correlation functions 
of W- with both w + and / are virtually flat lines, 
indicating the absence of any cross-correlations. 
The same is obtained for C[ui + , W^]. 

With this result in mind, consider again the 
partition function (26). We demonstrated before 
that we may essentially replace H R [W-,Q+] with 
i?[W_] and so rewrite Zqc as 



Z GC oc JvW-e- HGC JvQ + I GC [W-,u> + ]. 



(38) 



Since, as we have seen, / and W- are virtually 
uncorrelated, we argue that the results of the sim- 
ulation of W- cannot possibly be altered by the 
reweighting factor in a significant way. 
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FIG. 6. Correlation functions C[I,I] (dashed), 
C[(D+,d;+] (dashed-dotted), C[I,W-] (left inset), and 
C[Q+,W—] (right inset) in the same simulation as in 
Fig. 3. Note the different scale from Fig. 5. 

In sum, the results of these test simulations in- 
dicate that the W + fluctuations do not influence 
the structure of the polymer blends substantially. 
Therefore, the Monte Carlo simulations presented 
in the next section were restricted to W- fields 
only. 



IV. RESULTS AND DISCUSSION 

We have studied the two-dimensional ternary 
A+B+AB system at dimensionless number den- 
sity C — 50. In evaluating these simulations, we 
must establish two phase boundaries: first, the 
fluctuation-shifted order-disorder transition, and 
second, the onset of the phase-separated region. 
For the former, two parameters have been used 
which we shall define in the following. 

The Direction Persistence parameter. To define 
and determine this parameter, a given distribution 



W-(r) is translated via (15) into a "density" pat- 
tern 4>a(y)- All simulations performed with the 
field-theoretic Monte Carlo method that are pre- 
sented in this section took only W- fluctuations 
into account. Hence, 4>a{*) is real. Lattice points 
are denoted black if 0a (r) < 0.5, and white if 
4>a{y) > 0.5. We further define (max(Z a )) as the 
maximal length of either black or white sections 
along a one-dimensional cross section of the im- 
age in direction a, averaged over all offsets along 
the X axis (or equivalently the Y axis). The di- 
rection persistence parameter, A, is defined for a 
two-dimensional black-and-white image as 



A := 



( max G||)) 
(max(Zj_)} 



1, (39) 



all directions 



where l\\ denotes the direction being averaged over, 
and l± the direction perpendicular to a given /|| . A 
is in a disordered configuration and positive def- 
inite in a lamellar configuration, and thus consti- 
tutes a measure for the "lamellarness" of a config- 
uration. Furthermore, it is dimcnsionless and does 
not directly depend on the lattice size. A measures 
lamellarness from a local perspective. 

Anisotropy Parameters: As a complementary 
measure of anisotropy, we define the parameters 
F2 and F4,, which measure the anisotropy of the 
Fourier transform, F(q), of a 4>a{y) image to look 
at periodicity from a more macroscopic point of 
view: 

F n (q) ±-\ J%\F( q )\ 2 e^\. (40) 

By definition, Fi and F4 arc identically zero in a 
disordered (isotropic) configuration and are posi- 
tive otherwise. A more convenient way to analyze 
F n (q) is to calculate a normalized ratio of the in- 
tegral and the standard deviation of F n (q), which 
yields a single dimensionless number: 



F n := 



JdqF n (q) 



Jdqq 2 F n (q) ( JdqqF n (q) 
JdqF n (q) \JdqF n (q) 



r 



(41) 



1/2 



(42) 

F n , like F n (q), is for disordered and nonzero for 
lamellar configurations. Note that F(q) is not the 
structure factor of the melt but simply the Fourier 
transform of (j>A- In fully equilibrated simulations, 
the structure factor can be calculated from av- 
eraged expressions containing the fields W- and 
W+ AA . 




FIG. 7. Snapshots at C = 50, 4>h = 0, and 
= 1L1, 11.2, 11.3, 11.4, and 11.5. White points: 
< 4>a < 0.49, gray points: 0.49 < (j>A < 0.51, black 
points: 0.51 < 4>a < 1. All runs were started from 
disordered configurations. 



A. Monte Carlo simulations 

We shall first discuss the results of the Monte 
Carlo simulations. 

The shift of the order-disorder transition was ex- 
amined in the canonical ensemble on a 32 x 32 
lattice with dx=0.625 for 0# = and 0.2, and 
drr^O.78125 for <j) H = 0.4, 0.6, and 0.7. These 
dx values were chosen such that a single lamella 
was several pixels wide, which is a precondition for 
the evaluation of the A parameter to work prop- 
erly. For simplicity, we shall subsequently denote 
runs that started from a disordered configuration 
"D-started" and those that started out of a lamel- 
lar configuration, "L-started." In the latter, the 
lamellar periodicities of the starting configurations 
reflect the mean-field values at the given parame- 
ters. Moreover, in L-started runs, da; was chosen 
to allow for an integer number of lamellae to fit in 
the simulation box. 

In the case of copolymers only (<pH = 0), D- 
started runs at various \N yield the snapshots of 
Fig. 7, which were taken after approximately 1 mil- 
lion Monte Carlo steps and after the simulations 
were equilibrated. It is clearly visible from the im- 
ages that above a certain threshold %A^, a lamellar 
phase begins to form spontaneously. We did not 
observe any hysteresis effects. We should stress 
that the patterns 4>a (r) shown on these snapshots 
arc not to be confused with real density patterns 
- they are merely visualizations of the field distri- 
butions W-(r). 

To quantify the transition, we have plotted A 
in Fig. 8 and F2/F4 in Fig. 9. These plots also 
include curves for the other runs described above, 
i.e., homopolymer volume fractions 4>h up to 0.7. 
In all cases, a jump is displayed in both the A 
and F2 /F4 parameters at the fluctuation-corrected 
critical \N ■ For 4>h = 0, we find the transition 
to be shifted from the mean-field value of 10.495 
to 11.3(1), which is in good agreement with a re- 
sult obtained earlier by two of us with the complex 
Langevin method 43 . It is also further evidence that 



the W + fluctuations contribute only a minor cor- 
rection to the partition function. 

For homopolymer concentrations <j>n above 0.7 
the configurations did not spontaneously assem- 
ble into a lamellar phase for any %iV in D-started 
simulations on 32 x 32 lattices. The configurations 
can at best be described as very defective lamellae. 
However, in L-started runs the configurations only 
broke up at very low x^V values. For example, at 
4>h = 0.85, D-started runs were done for \N < 13, 
all of which stayed disordered, yet in L-started 
runs the lamellae only broke up for \N < 12.1. 
This obvious discrepancy can be made plausible by 
two arguments: First, in the D-started runs, as we 
approach higher \N values (which correspond to 
lower temperatures), we observe a freezing effect, 
i.e., the simulations get trapped in configurations 
with a certain degree of local periodicity yet the 
systems cannot move through configuration space 
efficiently enough to achieve quasi-long-range 
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FIG. 8. Direction persistence as measured by the 
A parameter vs. \N for different homopolymer vol- 
ume fractions <f>H at the order-disorder transition in 
simulations on 32 x 32 lattices. For 4>h = and 
0.2, the discretization da; = 0.625 was used, otherwise 
da; = 0.78125. 
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FIG. 9. F 2 (solid) and F 4 (dashed) vs. for dif- 
ferent homopolymer volume fractions (f>u at the or- 
der-disorder transition in the same simulations as in 
Fig. 8. 
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FIG. 10. F2 at homopolymer volume fraction 
<f>H = 0.7 (a) and 4>h = 0.8 (b) for different \N in sim- 
ulations on a 48 x 48 lattice that started from lamellar 
configurations. 

order. Second, in the L-started runs, in order to 
keep dx reasonably low (i.e., < 0.8), we can only 
put a small number of lamellae in the box. In the 
case of <pH = 0.85, that number is 3. But this 
in combination with the periodic boundary condi- 
tions imposed on the box in turn artificially stabi- 
lizes the lamellar phase, which is why we do not 
see a breaking up. Evidently, we must use bigger 
simulation boxes for <pn ^ 0.7. 

Because of the freezing effects at high \N de- 
scribed above, we restrict ourselves to L-started 
runs in examining the order-disorder transition on 
48 x 48 lattices. Even so, we find that the sim- 
ulations take an increasingly greater number of 
Monte Carlo steps to equilibrate the higher \N ■ 
For 4>h = 0.7, a time series of F 2 is displayed in 
Fig. 10 (a). The plot indicates that the transi- 
tion must be between %iV = 12.3 and 12.5: we 
see a substantial jump in both parameters, caused 
by a breaking-up of the initial lamellar configura- 
tion. This is in good agreement with the result 
from the 32 x 32 lattice. The F2 graphs at the 
transition, especially for %iV = 12.3, display a spe- 
cial feature: a periodic oscillatory motion once the 
run is equilibrated. The time scale of this oscil- 
lation can be explained as follows. Close to the 
transition, the lamellae do not easily dissolve com- 
pletely. Instead, defects develop which may open 
and close periodically, often alternating the con- 
tinuation of a given contour length between two 



defective lamellar branches. This oscillation is re- 
flected in F 2 . 



0.4 



For <f>H — 0.8, the corresponding time series is 
displayed in Fig. 10 (b). At x N = 12-5, the sys- 
tem is clearly disordered, and the transition is in- 
ferred to occur below \N = 13. Another run for 
(f>H = 0.85 indicates that that transition is at ap- 
proximately X-/V = 13, although even after 3 mil- 
lion Monte Carlo steps the configurations retained 
a strong correlation to the seed configuration. We 
conclude that on the one hand, using a bigger lat- 
tice has alleviated the artificial stabilizing effect 
of the periodic boundary conditions on the lamel- 
lae. On the other hand, however, we still encounter 
numerical freezing effects for \N ^ 13 at C=50, 
which prohibit the examination of this region of 
the phase diagram. 

As concerns the order of the ODT, we did sweeps 
in the grand canonical ensemble across the ODT 
for relative homopolymcr fugacitics z = 1 and z = 
2 (Fig. 11). Both plots are consistent with a very 
weakly first-order transition. In either case, the 
coexistence region has a width of no more than 
approximately 0.001 in 

The onset of the (macro) phase-separated two- 
phase region was examined in the grand canonical 
ensemble; a D-started run at parameters in the 
phase-separated region will eventually undergo 
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FIG. 11. Homopolymer concentration, vs. \N 
at z = 1 (a) and z = 2 (b) across the order-disorder 
transition. The dotted line indicates the xN coordi- 
nate of the ODT ±0.1. 
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FIG. 12. Difference of A and B monomer densities, 
A(j), vs. \N in a binary system of A and B homopoly- 
mers in the grand canonical ensemble. 
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FIG. 13. Difference in A and B monomer densities, 
A(j), vs. relative homopolymer fugacity, z, in the grand 
canonical ensemble. 



a spontaneous symmetry breaking and end up in 
a configuration dominated by either A or B seg- 
ments. Moreover, it is sufficient to do these runs 
on 32 x 32 lattices as no complex morphologies are 
investigated here. To quantify this behavior, we 
analyze the total difference of normalized A and B 
densities: 



A6 := 



7 



V- 1 / dr(<Mr) - <Mr)) 



(43) 



First, we looked at the demixing in a system con- 
sisting only of A and B homopolymers (4>h = 1)- 
In this case, Fig. 12 shows how the order pa- 
rameter A(j> grows from zero above approximately 
2.04/a, shifted up from the mean-field value of 
2/a. For various Fig. 13 displays A</> for 
relative homopolymer fugacities z (corresponding 
to (j)H < 1 in the canonical ensemble) around the 
transition. When translated into 4>h values in the 
canonical ensemble, the numerical accuracy ob- 
tained in the present work did not allow for de- 
termining the order of the corrected disorder/two- 
phasc transition. 
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FIG. 14. Anisotropy parameters obtained by Com- 
plex Langevin Simulations. The solid lines correspond 
to the F2 values while the dotted lines correspond to 
the F4 values. 
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FIG. 15. Averaged densities across the ODT, as ob- 
tained from Complex Langevin simulation runs. The 
homopolymer fraction is fixed at &h = 0.2. The 
XN values are: (a) x N = H- 4 ; (b) = U-7; ( c ) 
X N = 11.9. 



B. Complex Langevin simulations 

The same system has also been examined by 
complex Langevin simulations for selected param- 
eter values in the canonical ensemble. The sim- 
ulations were initiated from L-started configura- 
tions (in the lamellar phase) to examine the shift 
in the ODT. The transition points were identified 
by examining the values of the anisotropy param- 
eters F2 and F4. The behavior of these quantities 
(displayed in Fig. 14) is quite similar to that ob- 
served in the field-theoretic Monte Carlo approach, 
and exhibits a steep increase near the onset of mi- 
crophase separation. As an additional pictorial 
proof of such a transition, Fig. 15 displays the aver- 
aged values of the real component of the fluctuating 
density fields across the transition. 

As mentioned earlier in the text, complex 
Langevin simulations were hindered by numerical 
constraints arising from extensive phase fluctua- 
tions near the Lifshitz point, and required small 
time steps to maintain numerical stability. How- 
ever, the use of such time steps 

led to freezing effects, and rendered equilibration 
difficult. Consequently, we have displayed the 



results from CL simulations only for situations 
(i.e., <j)R values) for which the fluctuation-corrected 
transition temperatures could be identified reli- 
ably. It is encouraging to note that the results ob- 
tained from the CL simulations, the basis of which 
is somewhat distinct from the field-theoretic Monte 
Carlo method presented in this article, do indicate 
very good agreement with the results of the lat- 
ter approach. The implications of such results are 
threefold: (i) This serves to corroborate the results 
presented in earlier sections which were obtained 
using the Monte Carlo method, (ii) It also suggests 
that the fluctuation effects in the W + field, which 
were neglected in the Monte Carlo simulations (but 
which were included in the CL approach), have 
only a negligible influence upon the majority of 
the phase diagram, (iii) Further, the coincidence 
of CL results with the Monte Carlo approach sug- 
gests that despite the tenuous foundations of com- 
plex Langevin approach, such a method does yield 
physically reliable results for a major portion of 
the phase diagram, and thereby warrants further 
studies using such an approach. 



C. The fluctuation-corrected phase diagram 
and the formation of the microemulsion 

The results of our simulation investigations are 
put together in the fluctuation-corrected phase di- 
agram for C=50 (Fig. 16), and lead to a slope of 
the resulting fluctuation-corrected order-disorder 
transition which is in good qualitative agreement 
with the experimental result of Hillmyer et al. for 
a PEE/PDMS/PEE-PDMS system of intermedi- 
ate molecular weights 19 . In particular, the shape 
of the cusp-like region of bicontinuous microemul- 
sion identified therein is here reproduced compu- 
tationally for the first time. It is important to note 
that to minimize computational expense our sim- 
ulations were carried out in two dimensions, while 
the experiments arc inherently three-dimensional. 
In general, one would expect the fluctuation effects 
to be stronger in 2-d than in 3-d and, indeed, the 
microemulsion "cusp" appears broader in Fig. 16 
than in the experimental phase diagrams. Never- 
theless, we believe that a firm qualitative corre- 
spondence between simulation and experiment has 
been established. 
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preferential length scale of a configuration. Fig. 
17 compares Dc and Lq for various homopoly- 
mer concentrations 4>h at \N = 12.5, for both L- 
started and D-started runs. Whereas for lower 
Dc is larger than L , the two lengths cross each 
other as (pn is increased. Thus we have identified 
the mechanism by which fluctuations generate the 
microemulsion phase: it forms when the width of 
the lamellae grows larger than the curvature radius 
of the fluctuating boundaries, causing the lamellae 
to break up. This type of "Lindcmann criterion" 
for production of a bicontinuous microemulsion by 
melting of a lamellar phase has discussed elsewhere 
in the context of oil/water/surfactant phases 60 . 



FIG. 16. The fluctuation-corrected phase diagram 
at C=50. The circles show the results from the Monte 
Carlo simulations, the squares those from the complex 
Langevin simulations. The dotted lines indicate the 
mean-field result of Fig. 2. 
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FIG. 17. Curvature radius Dc (circles) and prefer- 
ential length scale, Lq, as obtained from the maximum 
of 5*o (?) (squares) at \N = 12.5 in simulations on a 
48 x 48 lattice. Filled and empty symbols correspond 
to disordered and lamellar initial conditions, respec- 
tively. Dashed lines are guides for the eye. 

To understand better what happens at the onset 
of the microemulsion phase, we have calculated the 
mean curvature diameter, Dc, of the boundaries 
of A and B microdomains: 



D c = 2 



Lcl 



ds 



n-l/2 



(44) 



where L c is the sum of all contour lengths of the mi- 
crodomain boundaries, and t is the tangent vector 
at a given coordinate s along the contour. Details 
of the calculation of this object will be presented 
elsewhere 59 . Furthermore, we examine the maxi- 
mum, q max , of Fo(q), which indicates the preferen- 
tial Fourier wave number of a configuration multi- 
plied by 2tt. In turn, L/q max =■ Lq indicates the 



V. SUMMARY AND CONCLUSIONS 



In the present paper, we have introduced a new 
way of evaluating the full partition function of a 
system of polymeric blends subject to thermal fluc- 
tuations. The complex Langevin approach which 
was employed in our earlier researches, accounts 
for both the W- and W+ fluctuations, and is able 
to predict the shift in the order-disorder transi- 
tion at low </>h and in the two-phase region. For 
0.7 < 4>h ^5 0.9 it runs into numerical difficul- 
ties. We believe that the large system sizes neces- 
sitated in such regions as well as the concomitant 
strong fluctuation effects arising near the Lifshitz 
points restricts the applicability of the complex 
Langevin method in such regimes. As an alter- 
native, we have proposed a field-theoretic Monte 
Carlo approach, which was shown to be applicable 
to greater regions of the phase diagram. Small lat- 
tices (32 x 32) sufficed up to <j)H — 0.7, above which 
bigger (48 x 48) lattices had to be used. In the re- 
gion where both methods work, they are in good 
agreement, indicating that the W + fluctuations ne- 
glected in the latter simulations contribute only a 
minor correction. Further evidence for this circum- 
stance was found in empirical test runs. Such ob- 
servations also suggest that one could have poten- 
tially employed a real Langevin scheme where the 
W+ fluctuations are neglected 50 . Regretfully, how- 
ever, a full analysis of the W+ fluctuations was ren- 
dered difficult by the occurrence of a non-positive- 
definite weighting factor (the sign problem). A 
fluctuation-corrected phase diagram was presented 
for C=50, which corresponds well in qualitative 
terms to the PEE/PDMS/PEE-PDMS system of 
intermediate polymerization indices examined in 
a series of recent experiments 16-19 , providing fur- 
ther corroboration of the existence of a region of 
bicontinuous microemulsion between the lamellar 
and phase-separated regions. The mechanism re- 
sponsible for the formation of this phase was shown 



to be intimately related to matching length scales 
of the curvature of microdomain boundaries and 
of the lamellar periodicity. In an upcoming arti- 
cle, we will examine properties of the disordered 
phase 59 . 
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APPENDIX: COMPLEX LANGEVIN 
METHOD 



In this appendix, we provide a brief outline of 
the complex Langevin method. Readers interested 
in more details are referred to Ref. 44 . The com- 
plex Langevin technique involves generating W- 
and W+ states in the entire complex plane (de- 
spite the fact that the integrations in Eq. (8) arc 
restricted to the real axis in W- and the imaginary 
axis in W+) by means of complex Langevin equa- 
tions. Specifically, the equations used to generate 
the real and imaginary parts of the W- and W+ 
fields, W R , Wi,W£, and W{ , are of the form: 

d[W R {r,t) + iWL(r,t)} _ SH[W-,W+] 



dt 



5W-(r,t) 
>(r,t), (45) 



d[Wl(r,t)-iW R (r,t)} _ ,5H[W-,W+ 



dt 



SW+(T,t) 

+0(r,i), (46) 



where 0(r, t) is a real Gaussian white noise with 
moments satisfying the appropriate fluctuation- 
dissipation relation: (9) = 0, (0(r, i)0(r', t')) = 
25(t — t')5(r — r'). In the above equations, 



SH[W-,W+] 
SW- (r) 

8H[W_, W+] 
SW+{r) 



C 



-4> A {v-w.,w+) 



2W- 
B (r;W-,W+) 



(47) 



C[4> A {v;W_,W + ) 

+ 4> B {v;W_,W + )-l] (48) 



46) produces in the long-time limit the generation 
of field states such that the "time" averages ac- 
quire the same values as those computed using the 
original non-positivc-defmite Boltzmann weight 48 . 
Equilibration of the simulation can be established 
by monitoring an appropriate quantity, such as 
the average of the imaginary component of en- 
ergy/densities which should identically vanish at 
equilibrium. We note that with such an averaging 
procedure, physical observables such as the aver- 
age A monomer density, <pA — (4>a), turn out to 
be real, while their fluctuating analogs such as 4>a 
are generally complex. 

It is pertinent to note a few special features of 
the complex Langevin equations (45 - 46). In the 
absence of the noise term 0(r, t), the dynamics will 
drive the system to the nearest saddle point. Con- 
sequently, it is straightforward to probe the sta- 
bility of the mean-field solutions within the CL 
simulation scheme by setting the noise term identi- 
cally to zero. Furthermore, this feature enables the 
sampling in the presence of the noise to be auto- 
matically confined near a saddle point, avoiding a 
long equilibration period if the system were to be 
prepared in an initial configuration far from the 
saddle point. 



The evolution of the W R , WL and Wf , W{ fields 
according to the above Langevin equations (45 - 
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